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I. Introduction 

A DJOINT solutions of the governing flow equations are becoming increasingly important for the devel- 
opment of efficient analysis and optimization algorithms. A well-known use of the adjoint method 
is gradient-based shape optimization. 1,2,3,4 ' 3 ’ 6 Given an objective function that defines some measure of 
performance, such as the lift and drag functionals, its gradient is computed at a cost that is essentially 
independent of the number of design variables (geometric parameters that control the shape). More recently, 
emerging adjoint applications focus on the analysis problem, where the adjoint solution is used to drive mesh 
adaptation, ' ,8,9, 10 as well as to provide estimates of functional error bounds and corrections. 11 The attrac- 
tive feature of this approach is that the mesh-adaptation procedure targets a specific functional, thereby 
localizing the mesh refinement and reducing computational cost. 

Our focus is on the development of adjoint-based optimization techniques for a Cartesian method with 
embedded boundaries. 12 In contrast to implementations on structured and unstructured grids, Cartesian 
methods decouple the surface discretization from the volume mesh. 13 This feature makes Cartesian meth- 
ods well suited for the automated analysis of complex geometry problems, and consequently a promising 
approach to aerodynamic optimization. Melvin et al. 14 developed an adjoint formulation for the TRAN AIR 
code, 15 which is based on the full-potential equation with viscous corrections. More recently, Dadone and 
Grossman 16 presented an adjoint formulation for the Euler equations. In both approaches, a boundary con- 
dition is introduced to approximate the effects of the evolving surface shape that results in accurate gradient 
computation. 

Central to automated shape optimization algorithms is the issue of geometry modeling and control. The 
need to optimize complex, “real-life” geometry provides a strong incentive for the use of parametric-CAD 
systems within the optimization procedure. In previous work, 17 we presented an effective optimization 
framework that incorporates a direct-CAD interface. 18 In this work, we enhance the capabilities of this 
framework with efficient gradient computations using the discrete adjoint method. We present details of the 
adjoint numerical implementation, which reuses the domain decomposition, multigrid, and time-marching 
schemes of the flow solver. Furthermore, we explain and demonstrate the use of CAD in conjunction with the 
Cartesian adjoint approach. The final paper will contain a number of complex geometry, industrially relevant 
examples with many design variables to demonstrate the effectiveness of the adjoint method on Cartesian 
meshes. 
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II. Problem Formulation 


The aerodynamic shape optimization problem consists of 
determining values of design variables X. such that the ob- 
jective function J is minimized 


min J(X,Q ) 

A 

subject to constraint equations Cj 

Cj(X,Q)< 0 j = 1, . . . , N c 


(1) 


( 2 ) 
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where the vector Q denotes the continuous, conservative flow 
variables and N c denotes the number of constraint equations. 
The flow variables are forced to satisfy the governing flow 
equations within a feasible region of the design space P 
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Figure 1. Components of the analysis module 


which implicitly defines Q = f(X). 

The objective function defines the goals of the optimiza- 
tion problem, while the constraint equations limit the feasi- 
ble region of the design space. The constraints may involve 
performance functionals, such as lift, geometric quantities, 
such as volumes and thicknesses, and also simple bound con- 
straints for design variables. Performance objectives can be 
specified for the entire configuration or for a specific sub- 
set of components, for example moments on control surfaces. 

This is accomplished by' using the Geometry Manipulation 
Protocol, 19 where we specify component hierarchies for parametric-CAD assemblies and triangulations to 
intuitively reflect aerodynamic design goals. 

A modular framework 17 is used to solve the optimization problem defined by Eqs. 1-3. We cast the 
optimization problem as an unconstrained problem by lifting the side constraints, Eq. 2, into the objective 
function using a penalty' method. The constraint imposed by the flowfield equations, Eq. 3, is satisfied at 
every' point within the feasible design space, and consequently these equations do not explicitly' appear in the 
formulation of the optimization problem. An unconstrained BFGS quasi-Newton algorithm coupled with a 
backtracking line search 20 ' 21 is used to find the optimal solution. At the core of the optimization framework 
is the analysis module, which consists of a CAD sy'stem interface provided by CAPRI, 22,18 a Cartesian grid 
generator for component-based geometry, 23 and a flow solver, as outlined in Fig. 1. Below, we provide a 
brief description of the flow solver in order to clarify the subsequent discussion. Thereafter, we focus on the 
development of the gradient computation algorithm. 


III. Governing Flow Equations and Numerical Method 

The governing flow equations are the three-dimensional Euler equations of a perfect gas. For a finite 
region of space with volume V and surface area S, the integral form of the Euler equations is given by 


d t 


Q dV 


F ■ n dS — 0 


( 4 ) 


2 of 12 


American Institute of Aeronautics and Astronautics 




k 


(5) 


where Q = [p, pu. pv. pw, pE] T , F is the inviscid flux tensor 
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and n is the outward facing unit normal vector. 

The Euler equations are solved with a finite- volume method on a regular Cartesian mesh with embedded 
boundaries. The mesh consists of hexahedral cells, except for a layer of body-intersecting cells, or cut-cells, 
that are arbitrary polyhedra adjacent to the boundaries. 23 Spatial discretization is based on a cell-centered 
approach, where the control volumes V correspond to the mesh cells and the cell-averaged value of Q, denoted 
by Q, is located at the centroid of each cell. The control volumes are fixed in time, resulting in the following 
semi-discrete form of Eq. 4: 

/>§ + «(<?)=» ( 6 ) 

where Q = [Qi> Qi, ■ . . , <5.v] T is the discrete solution vector for all N cells, D is a diagonal matrix containing 
the corresponding cell volumes, and R is the residual vector. The residual in each cell i is expressed as 

R i = J2 F i' n i S i ( 7 ) 

j£Vi 

■where j denotes the jth. face of volume V) with area S. The flux is evaluated at the face centroids based on 
a midpoint quadrature rule. 

The residual evaluation for a second-order accurate discretization proceeds by first reconstructing the 
solution to the cell face, given by 

Ql = Qi + d L 4>iX?Qi ( 8 ) 

where Ql denotes the left state at a given face reconstructed from cell i, is the distance from the cell 
centroid to the face centroid, d> is the limiter (in slope form) used to enforce monotonic solutions, and V<3 is 
the solution gradient determined via a linear least-squares procedure. A similar formula is used to reconstruct 
the right state, denoted by Q r. Next, the flux-vector splitting approach of van Leer is used to resolve the 
flux discontinuity at the face, represented by 

F(<2l,<3r) ■ n = f + (Q L ,™) + f “(QR, n ) (9) 

The assembly of the residual vector is accomplished by a loop over the faces of the mesh, where the flux 
contributions axe scattered from the faces and accumulated in the cells. Steady-state solutions are obtained 
using a five-stage Runge-Kutta scheme in conjunction -with multigrid and domain decomposition for parallel 
computing. For further details on the flow solution algorithm, see Aftosmis et al } 2 - 24 and Berger et al. 2b 

IV. Adjoints and Sensitivities 


A. Formulation 

The gradient, G. of the discrete objective function J \X.Q(X) 


bv 


dj = dj_ 5JdQ 
dA ~ 3X + qq dX 
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where Q is the steady-state solution of Eq. 4 for a given set of design variables 

R{X,Q) = 0 (11) 

Note that this is the discrete equivalent of Eq. 3. We reduce the vector of design variables, X, to a scalar 
in order to clearly distinguish between partial and total derivatives. For problems with multiple design 
variables, it may be helpful to note that G and dJ/dX are [1 x A'd] row vectors, dJ/dQ is a [1 x Np] 
row vector, and dQ/dX is a [Np x A r n] matrix, where Np, and Np represent the number of design and flow 
variables, respectively. 

In Eq. 10, the evaluation of the term dQ/dX, referred to as the flow sensitivities, is obtained by differ- 
entiating Eq. 11 with respect to the design variables 

dX dX dQdX 1 ’ 


where we assume that the implicit function Q(X) is sufficiently smooth. 26 ’ 4 5 Realizing that = 0, since 
for any design variable Eq. 11 is always satisfied, Eq. 12 simplifies to the following large sparse system of 
linear equations 

dRdQ__dR 

dQdX dX { ’ 

The direct, or flow-sensitivity, method results from solving Eq. 13 for the flow sensitivities dQ/dX and using 
these values in Eq. 10 to obtain the gradient. 

In order to formulate the discrete-adjoint method, substitute Eq. 13 into Eq. 10 to obtain 


dj _ dj dj_ fdR\ 1 dR 
d X~dX~d$\dQJ dX 


From the triple-product term in Eq. 14, define the following intermediate problem 

BQ \dQ J 


(15) 


where ip is a [Np x 1] column vector. Post-multiplication of both sides by dR/dQ and applying the transpose 
operator results in the following linear system of equations 


dR T -._ 
dQ dQ 


(16) 


This is known as the adjoint equation, and the vector ip represents the adjoint variables. Substituting ip into 
Eq. 14, the expression for the gradient becomes 


d J _ dj - T dR 
dX ~ dX ~ ' dX 


(17) 


Note that Eq. 13 is a linear sj^stem with multiple right-hand sides dependent on the number of design 
variables. In contrast, Eq. 16 is independent of the design variables, which is the reason for the efficiency of 
the adjoint method. 
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B. Numerical Implementation 


Our goal is to reuse the solution method of the flow solver, in particular the Runge-Kutta time marching, 
local time stepping, multigrid, and domain decomposition schemes, to solve the adjoint and flow-sensitivity 
equations. This minimizes the code development time and simplifies code maintenance issues. A promising 
approach for generating flow-sensitivity (linear) and adjoint codes is automatic differentiation. 2, ,2S ’ 29,30 
Although we do not use automatic differentiation in this work, we found that understanding the underlying 
principles and techniques is helpful for producing accurate and efficient sensitivity codes. Overall, the 
resulting approach is similar to the work of Giles et al. 31 

The solution of the adjoint and flow-sensitivity equations proceeds by introducing unsteady terms in 
Eqs. 13 and 16 to obtain 


D** + ™Q> + ™ =0 

d* [dQ^ dX] 

d < \d Q dQ J 


(18) 

(19) 


■where Q’ = dQ/dX. A key step in the solution procedure is the evaluation of the matrix-vector product 
terms 

and d -£-i (20) 

dQ dQ 

We consider the flow-sensitivity equation first, since the explicit linearization of the residual equations can 
be avoided. The matrix-vector product is approximated as follows 


dR~, R (Q + eQ') - R(Q) 

— — 

dQ e 


( 21 ) 


with e = a/ (e m )/||v|j. 32,6 The final paper will include a detailed discussion of the flow-sensitivity and adjoint 
solvers. 

The partial derivative term dj /dQ in Eqs. 10 and 16 is evaluated using finite differences. The remaining 
partial derivative terms in Eqs. 10 and 17, namely the objective function sensitivity dj fdX and the residual 
sensitivity dR/dX, are usually approximated with finite differences as well. Although in some instances the 
partial derivatives can be obtained analytically, finite differences are necessary -when a CAD system is used 
to parameterize the geometry. For example, the centered-difference formula for the residual sensitivities is 
given by 


gR R(X + he n ,Q ) 

dX n = 


R X - he r 


,Q) 


2 h 


( 22 ) 


where e n denotes the nth unit vector, 


h = max (eo • \X n \, (23) 

and n — 1, ... ,Nu- Typical value of eo is 1 x 10 -3 . It is important to realize that Eq. 22 involves two 
evaluations of only the residual vector per design variable and not two flow solutions, i.e. Q is constant. 

Before presenting our numerical implementation for the computation of dj /dX and dR/dX on Cartesian 
grids, it is insightful to briefly review the implementations on structured and unstructured grids. The 
standard finite-difference approach relies on a grid-perturbation strategy to smoothly “convect” the nodes of 
the volume grid w'ith the perturbed surface geometry, as sketched in Fig. 2. This results in an accurate and 
efficient evaluation of the partial derivative terms. .Alonso et al. 33 have successfully extended this approach 
to a CAD-based design environment by introducing geometry patches wuthin the parametric-CAD model. 
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An important observation is that the grid-perturbation strat- 
egy explicitly couples shape sensitivities to volume grid sensitiv- 
ities. However, a grid-perturbation without a corresponding sur- 
face perturbation should not, in principle, influence the objective 
function gradient. On the basis of this argument, Jameson and 
Kim 34 introduced a reduced adjoint gradient formulation where 
only surface perturbations are considered. Their results indicate 
that for shape optimization problems the reduced approach works 
well. Similar results are also presented by Soto and Lohner. 35 In 
contrast, work by Anderson and Venkatakrishnan 3 and later by 
Nielsen and Anderson 36 shows that for cases that involve rigid 
body motion, for example flap position optimization, the use of 
a grid-perturbation strategy is important. This appears to be re- 
lated to the treatment of the trailing-edge singularity. 

Turning our attention back to embedded-boundary Cartesian 
grids, we note that a perturbation of the surface shape affects 
only a few near-by cells of the volume grid, see Fig. 3. Therefore, 
the extent of grid sensitivities is naturally limited to roughly 0(N 2 ) cells for a volume grid with (D(N 3 ) 
cells. However, a straightforward implementation of finite differences is complicated by the emergence and 
disappearance of cut-cells at the wall boundary, as well as the relative motion of the discretized surface with 
respect to the volume grid. 

To circumvent these difficulties, we use the following approach. 

For dJ/dX, the converged flow solution is reconstructed to the 
vertices of the underlying surface triangulation. Given a small 
perturbation, a new surface triangulation is generated using the 
CAPRI interface. We force the new triangulation to have a one- 
to-one triangle mapping with respect to the baseline triangulation. 

This is accomplished by using the parametric values of the base- 
line triangulation when tessellating the perturbed surface. The 
objective function values required for the finite-difference approxi- 
mation are now evaluated using the surface triangulations instead 
of the volume grid. Hence, the potentially non-smooth changes in 
the cut-cells are avoided, and the enforcement of constant Q in the 
evaluation of this partial derivative term becomes a trivial task. 

This approach works well for both shape and rigid-body motion 
design variables. 

For the remaining term dR/dX, we have implemented and 
studied two approaches, both based on finite-differences. In the 
first approach, we assume that the only parameter dependent on the design variables in Eq. 11 is the surface 
normal. This approach is similar to the “transpiration” boundary condition used in TRANAIR. 37, 14 It is 
relatively easy to implement, and the evaluation of dR/dX is fast since the volume grid remains unchanged. 
The accuracy of this approximation can be sufficient when the design variables involve shape changes normal 
to the surface including wing-twist distributions. However, a possible limitation of this approach may be 
rigid-body motion. 

In the second approach, we construct a new volume grid for the perturbed surface triangulation. The 
baseline solution is interpolated to the new grid via a fast interpolation algorithm presented in Ref. 24. We 
evaluate the residual on the new grid, which is then passed back to the baseline grid for the finite-difference 
computation. Cut-cells that appear or disappear during this procedure are tagged, and we either zero their 
contribution to the residual sensitivity computation, or we interpolate a value from their face neighbors. We 
evaluate both approaches in the following section. 



Figure 3. Cartesian grid with a baseline 
(red) and perturbed (green) airfoils 



Figure 2. Grid-perturbation approach 
(baseline airfoil and grid are red, pertur- 
bation is green) 
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V. Results and Discussion 


In this abstract we study a simple, two- 
dimensional design example to validate our ap- 
proach. The baseline geometry is the NACA 0012 
airfoil parameterized with a cubic B-spline curve 
using 17 control points. We choose the vertical 
motion of a control point near the leading edge of 
the airfoil as the design variable, shown in Fig. 4. 

The freestream Mach number is 0.5 and the an- 
gle of incidence is fixed at 2 deg. The goal of 
the optimization is lift enhancement. The base- 
line Cl is 0.245 and we specify a target Cl of 
0.5. The “volume” grid contains roughly 9,000 

cells, but since this is a two-dimensional problem, , _ 

, , . . , . . , Figure 4. Design variable and B-spline parameterization 

the depth m the span-wise direction contains only of the NACA 0012 airfoil 

two cells. To help validate the results, we solve a 

similar problem using Optima2D , 6 which is an aerodynamic optimization tool based on a structured-grid 
approach. 

We begin by examining the computation of the dJ/dX term in Eq. 10. We evaluate the partial derivative 
of the lift functional (objective function) for the above problem with the angle of incidence set to 0 and 2 
deg. Note that at 0 deg., this partial derivative term should vanish, providing a good validation exercise. 
The result at 0 deg. is indeed zero, and at 2 deg. we obtain —2.81 x 10~ 3 . Optima2D predicts —1.4 x 10“ 2 , 
which is a reasonable agreement considering the differences of the two approaches. 

Prior to solving the flow- sensitivity equation, we evaluate the residual sensitivity dR/dX. Figure 5(a) 
shows dR/dX values for the y-momentum equation. These are obtained by constructing (or re-cutting) a 
new volume grid for the perturbed shape and interpolating the steady-state solution. Figure 5(b) shows the 
results when only the surface normals are changed. The agreement in Fig.5 is excellent; however, we note 
that for the remaining field variables significant errors may arise due to this approximation. We plan to 
investigate this further in the final paper. 

The solution of the flow-sensitivity equation is shown in Fig. 6. Note that the zone of influence for 
the shape design variable extends from the leading edge to roughly 60% chord on the upper surface. The 
convergence of the flow and sensitivity equations is shown in Fig. 7. The flow equations converge in 500 
multigrid cycles (only a 2-level multigrid is used), and the flow sensitivities require additional 400 cycles. 
The final gradient value is —1.35. Optima2D computes a value of —1.67, while a finite-difference estimation 
of the gradient by recomputing the flow gives —1.55. A grid-refinement study will be included in the final 
paper to validate gradient accuracy. Furthermore, the final paper will contain additional design examples to 
demonstrate the effectiveness of the proposed approach. 
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(a) Sensitivity of Eq. 11 via the grid re-cutting and interpolation approach 



(b) Sensitivity of Eq. 11 using variations of only surface normals 


Figure 5. Contours of 6R/8X for the y-momentum equation 


8 of 12 



American Institute of Aeronautics and Astronautics 








(b) Structured-grid approach using OptimaSD 6 

Figure 6. Contour plots of flow sensitivity: changes in density with respect to a shape design variable 
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Figure 7. Convergence to steady-state for the mean flow and the flow-sensitivity equation 
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